Printed 3 January 2013 



Molecular dissipation in the nonlinear eddy viscosity in the 
Navier-Stokes equations: modelling of accretion discs 

G. Lanzafame* 

INAF - Osservatorio Astrofisico di Catania, Via S. Sofia 78 - 95123 Catania, Italy 



Accepted -. Received ; in original form 

ABSTRACT 

Physical damping, regarding the nonlinear Navier-Stokes viscous flow dynamics, refers 
■ to a tensorial turbulent dissipation term, attributed to adjacent moving macroscopic 

flow components. Mutual dissipation among these parts of fluid is described by a 
braking term in the momentum equation together with a heating term in the energy 
equation, both responsible of the damping of the momentum variation and of the 
viscous conversion of mechanical energy into heat. 

A macroscopic mixing scale length is currently the only characteristic length 
" 53 ■ needed in the nonlinear modelling of viscous fluid dynamics describing the nonlin- 

ear eddy viscosity through the kinematic viscosity coefficient in the viscous stress 
tIh ■ tensor, without any reference to the chemical composition and to the atomic dimen- 

sions. Therefore, in this paper, we write a new formulation for the kinematic viscosity 
coefficient to the turbulent viscous physical dissipation in the Navier-Stokes equations, 
where molecular parameters are also included. 

Results of 2D tests are shown, where comparisons among flow structures are made 
on 2D shockless radial viscous transport and on 2D damping of collisional chaotic 
turbulence. An application to the 3D accretion disc modelling in low mass cataclysmic 
variables is also discussed. 

Consequences of the kinematic viscosity coefficient reformulation in a more strictly 
physical terms on the thermal conductivity coefficient for dilute gases are also dis- 
cussed. 

The physical nature of the discussion here reported excludes any dependence by 
the pure mathematical aspect of the numerical modelling. 

Key words: accretion, accretion discs - hydrodynamics - binaries: close - stars: 
novae, cataclysmic variables. 



1 INTRODUCTION 

Physical dissipation in the viscous fluid dynamics is the only 
physical mechanism for the attenuation of the momentum 
transfer and for the conversion of mechanical energy (kinetic 
+ potential) into heat. It originates from microscopic parti- 
cle interactions on molecular scale lengths. For this reason, 
it is not currently included in the nonlinear Navier-Stokes 
equations for viscous flows, where macroscopic spatial reso- 
lution lengths of moving fluid components are much larger 
than molecular scale lengths. 

Thus, a physical turbulent viscosity is used in the 
Navier-Stokes equations, as a tensorial viscous dissipation 
term, relative to mutual interactions among contiguous 
macroscopic moving flow parts, producing their braking and 
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their contemporary heating. A turbulent kinematic viscos- 
ity coefficient is characterized by a macroscopic scale length 
multiplied by a scale velocity in the Von Karman descrip- 
tion, originally formulated to describe a repeating pattern 
of swirling vortices caused by the unsteady separation of 
flow of a fluid over bluff bodies. Hence, a mixing length is 
often required in the formulation of the kinematic viscos- 
ity coefficient v in the viscous stress tensor describing the 
nonlinear turbulent eddy viscosity, without any reference to 
the chemical composition and to the microscopic molecular 
dimensions. Moreover, physical turbulent viscosity often in- 
cludes arbitrary para meters, to be se t case by case, as it is 
for of the well known IShakural (| 19721 ); IShakura fc Sunvaevl 
formulation for disc structures. 

A dissipation mechanism is always necessary in the com- 
putational collisional fluid dynamics, even in the non viscous 
modelling to solve the strictly hyperbolic Euler equations, if 
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flow discontinuities (the Riemann problem) must be solved. 
In the physically inviscid fluid dynamics, "shock captur- 
ing" methods adopt either an artificial viscosity contribu- 
tion or take advantage of some moderate intrinsic numeri- 
cal dissipation. Instead, "shock tracking" methods develop 
a dissipation separately handling shock fronts using appro- 
priate Riemann solver algorithms, through algebraic aver- 
ages between the two left-right sides as Rankine-Hugoniot 
jump conditions, in the Godunov-type methods. In addi- 
tion, a further dissipation is unavoidably in t rinsically gen - 
erated by truncation errors (jFletcher! [l99ll ; Iffirschl ll99Tt) . 
In the finite difference methods, dissipation comes about 
by second order derivatives coming from the Taylor series 
expansion f or incremental ratio s of the first order spatial 
derivatives ()Park fe Kwonll2003h . especially for implicit in- 
tegration techniques. Such a dissipation contribution is also 
useful to smooth out spurious heating and to treat transport 
phenomena. 



The physical, turbulent dissipation and the fictitious ar- 
tificial or numerical ones are conceptually di s tinct, although 
forma ll y similar, as shown inlMolteni et"aD (1991); Murray! 
1996); Okazaki et al.l (|2002^ . and discussed by Lanzafama 



2008 . 12009) in the case of smooth particle hydrodynamics 



(SPH) accre tion disc modelling in close binaries (CBs), or by 
iTord [| 19991 ) in the case of finite difference or for finite vol- 
ume integration techniques. In both cases, some arbitrary 
parameters, to be tuned case by case, are included and/or a 
dependence on the spatial resolution length affects the nu- 
merical dissipation. 

Therefore, in this paper we propose the formulation of 
a physical turbulent kinematic viscosity coefficient in the 
Navier-Stokes equations, free of any arbitrary parameters, 
where microscopic physical characteristics are also included. 
A macroscopic scale length (mixing length) is obviously still 
used, since solutions of the nonlinear Navier-Stokes equa- 
tions, involve macroscopic physical properties. At the same 
time, a reformulation for v in a more physical sense is also 
coherent to a reformulation for the thermal conductivity co- 
efficient c for dilute gases. 

In §2 of this paper we discuss some general aspects of 
dissipation in the computational collisional fluid dynamics; 
in §3, we shortly describe characteristics of some adopted 
physical turbulent kinematic viscosity coefficients v, while 
in §4 we formulate the physical development of both v and c 
coefficients including microscopic molecular characters. In §5 
we show some results for some essential 2D tests on shockless 
radial viscous transport in an annular ring, as well as for the 
damping of 2D Burger's turbulence. Instead, in §6, we show 
an astrophysical application in the case of a 3D accretion 
disc modelling in a low mass close binary (LMCB) system, 
comparing 3D accretion disc structures obtained by using 
different v coefficients and gas compressibility. 

Despite the adoption of a numerical technique intrinsi- 
cally quite viscous, like the SPH, successful physically vis- 
cous results unconditionally show that the new formulation 
for the kinematic viscosity coefficient works well without any 
restriction. 



2 DISSIPATION IN VISCOUS AND NON 
VISCOUS FLUID DYNAMICS 

In the physically non viscous flows, the hyperbolic Eulcr 
system of equations 

dp . . . , 

— + p\/ ■ v = continuity equation(l) 

dv Vp „ . 

— = h J momentum equation(z) 

dt p 

-7- ( e + -v 2 \ = V ■ (pv) + f ■ v energy equation(3) 

dt \ 2 / p 



dr 



kinematic equation. (4) 



must be solved, together with the state equation (EoS) 
of the fluid 



state equation(5) 



Most of the adopted symbols have the usual meaning: 
d/dt stands for the Lagrangian derivative, p is the gas den- 
sity, e is the thermal energy per unit mass, p is the ideal 
gas pressure, here generally expressed as a function of local 
properties, v and r are the vector velocity and position, / is 
the external force field per unit mass. The adiabatic index 7 
has the meaning of a numerical parameter whose value lies 
in the range between 1 and 5/3. 

Since the Riemann problem must be correctly solved for 
collisional flows in the case of shocks, a dissipation mecha- 
nism is necessary otherwise frontal colliding flows trespass 
each other. Such a dissipation could be either explicit, as 
an artificial viscosity term for shock capturing schemes, es- 
pecially for finite volume schemes, or it c ould be intrin- 
sic through a s pecific Riemann solver code jLeVeaueill992l; 
iFletcherl Il99ll ; IffirsdJ Il997l ; iLeVeaud |2002| ; iPark fc Kwonl 



120031 ) either for shock tracking schemes or for Eulerian finite 
difference schemes by commuting mathematical derivatives 
in incremental ratios. As an example, in the finite difference 
techniques, the conversion of the 1st order spatial derivative 
of the generic physical quantity u is (m — 114-1)/ Ax, where 
i is a spatial grid index and Aa; is the grid spatial resolution 
length. For a better stable result, the same incremental ratio 
is rewritten as 

m — Uj-i _ Uj+i — Uj-i Uj+i — 2uj + Uj-i . . 

Ax ~~ 2A~i 2Ai ' ( ' 

which gives a much better stability, at the cost of a 
reduction in accuracy. The contribution of the second term 
(ui+i — 2ui +Ui-i)/ Ax analytically corresponds to a 2nd or- 
der spatial derivative, working exactly like a real viscous non 
physical contribution. These manipulations of spatial deriva- 
tives in the incremental ratios in finite terms are necessary 
to ensure stability to the solutions of hyperbolic systems of 
equations. This means that the inclusion, or the numerical 
development of a non physical dissipation distorts numerical 
results especially for non collisional events like shear flows 
or transport phenomena. These numerical difficulties arise 
when the EoS 

p — (7 — l)pe perfect gas equation(7) 

is adopted for ideal flows. Instead, an EoS as: 



-1/3 



V • v 



3c 



(8) 
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includes a real macroscopic physical dissipation cor- 
rectly handling the Riemann problem, as well a s transport 
and sh ea r flows free of any local gas compression (|Lanzafamd 
l2010al lbl; [Lanzafame et~afll201ll '). c s is the sound velocity, n 
is the numerical density, while 



C 



- arccot 



(9) 



where D>1 and where vr is the component of velocity 
along the direction of collision. D is a large number describ- 
ing how much the flow description corresponds to that of an 
ideal gas: D ~ A/d, being A oc p~ 1//3 the molecular mean 
free path, and being d the mean linear dimension of gas 
molecules. The physical dissipation, expressed by the two 
further terms in eq. (8) (the linear and the quadratic terms 
in V • v) of the reformulated EoS, better treats both shocks 
and shear flows, even in a Lagrangian description. Their 
inclusion substitutes artificial viscosity terms and does not 
represent a physical turbulent viscous contribution, but the 
real physical dissipation coming out because eq. (7) of the 
EoS should strictly be applied only to macroscopic static 
or quasi-static processes. Notice that in this case, this real 
macroscopic physical dissipation does not originate from a 
physical viscosity. Instead, it originates from the irreversible 
thermodynamic process and is better evidenced in a La- 
grangian description. 

In the physically viscous flows, the Navier-Stokes equa- 
tions explicitly include macroscopic physical dissipation 
terms in the momentum and in the energy equations: 



dv 
~dt 



p p 



3 THE PHYSICAL KINEMATIC VISCOSITY 
COEFFICIENT 

Typical kinematic laboratory viscosities are of the order of 
v = 0.001 - 1 cm 2 s" ' 
in the ratio 



to be compared with inertial forces 



Re = 



inertial forces It 



flow^flow 



viscous forces 



(14) 



where lfi ow and Vfi ow are the characteristic length and 
velocity scales of the microscopic flow. Laboratory experi- 
ence shows that for Re > Re cr it ~ 10 2 — 10 3 , flow becomes 
turbulent. Re cr u is the characteristic Reynolds number as 
observed so far. 

In the full nonlinear approach, the full non linearity 
of the Navier-Stokes equations is considered, where spatial 
derivatives of the entire velocity field are used. Neither the 
Reynolds averages of the Navier-Stokes equations in boxes 
of intermediate size (as in the linear approach), nor the full 
Navier-Stokes equations working with spatial gradients of 
the me an velocity fiel d (as in the nonlinear Boussinesq ap- 
proach (|Schmittl(200"7h ), are considered. 

To characterize a nonlinear macroscopic physical kine- 
matic viscosity coefficient v, characteristic length and veloc- 
ity scales I and v are needed, which are unknown, in princi- 
ple. 

Typically (|Prandtllll925h , a mixing length model can be 
used, where 



dv 

d-x 



dv 

dx 



(15) 



(16) 



Navier-Stokes momentum equatio^lO) 



or, more generally, for a better statistical evaluation, 



dt V + 2 J = ~p V ' ^ PV ~ V ' T) + cV(pe)] + f ' V 

Navier-Stokes energy equation, (11) 

where the viscous stress tensor r and the thermal con- 
ductivity 1/pV ■ [cV(pe)] terms are explicitly added, to be 
solved together with the continuity equation, the kinematic 
equation and the EoS. It is important to note that the ther- 
mal flux term includes two contribution: the first contribu- 
tion depends on the thermal gradient (Ve) that is currently 
used for solids or for incompressible fluids, while the second 
one depends on the density gradient (Vp) correlated to the 
mass diffusivity, here not included in the continuity equa- 
tion. The matrix element (a, /3) of the viscous stress tensor 



T a ,p = V a &,P + CV ■ u 
and 

dv a dvp 2 

Oq,/3 = h » -Oa^V • V. 

OX/3 OX a 3 



(12) 



(13) 



r) and C, are the dynamic first (shear) and second (bulk) 
physical viscosity coefficients. 

In the present study, we simply consider C, = and eq. 
(7) as EoS. By definition, the physical kinematic viscosity 
coefficient is v = rj/p. 



dv x 
dy 



dv x 
dz 

dv z dv z 
dx dy 



+ 



I dv v 



8Vy^ 2 



\ dx 



+ 



(17) 



Here, the problem relies in the evaluation of I. Being 
h the computational spatial resolution, and being L the 
scale length of the entire computational domain, h ^ I ^ 
L. Because of the lack of any geometric information, the 
only physical scale lengths we know are those relative to 
the hydrostatic equilibrium (in the presence of an external 
force field): J dp/pf, as well as p/\\7p\, p/|Vp|, |w|/V • v, 
|u|/|V x v\, etc.. 

Since, in 3D the natural tendency is the de velopment of 
smalle r structures in a d irect cascade process feolmogorovl 
Il941al lbh. some authors (|Trampedach fc Steinlfeoilj ) calcu- 
late 



(18) 



where U refers to various scale lengths as (d\n p/dr) 1 , 
(dlnV/dr)- 1 , (dlnp/dr)-\ (V ■ v/\v\)-\ etc.. 

Instead, no information we have about v, since we only 
know Vfi ow and c s . 

In the vi scous acc r etion d isc modelling, the Shakura 
and Sunyaev IShakural (|l972f ); IShakura fc Sunvaevl (| 19731 ) 
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parametrization of turbulent viscosity is largely adopted. In 
this approach, the kinematic viscosity coefficient is 



-Iv, 



(19) 



where both I and v are unknown. Assuming the flow 
isotropy, a Keplerian tangential kinematics, and the vertical 
hydrostatic equilibrium, 



I = a t H, 



(20) 



where H ~ rc s /vKe P i, the local disc thickness, is the 
local shortest macroscopic scale length, and Q; ^ 1 is a 
scaling quantity. Without any isotropy assumption, ai > 1. 
At the same time, 



v = a v c s 



(21) 



where a v ^ 1. Whenever v > c 3 , shocks would dissipate 
the energy, reducing the velocity to subsonic. Hence, in the 
Shakura and Sunyaev approach, 



v = ^ona v c a H = assCsH, 



(22) 



with ass < 1 to be found. 

IPringld |l98lft found 0.01 sC a SS 0.03 for 0.1 ^ 
£/£ E ^ 1 (jC-e = Eddington luminosity) as a lower limit for 
active galactic nuclei (AGN). For numeri cal simulations of 
AGN , ass ~ 10 -4 — 10~ 3 is often adopted (|Lanzafame et ail 
1 19981 . 120081 1. Values for a S s « 10~ 2 have also been foun d 
for the observed protostellar objects (|Hartman et al.ll 1998h . 
while for a fit of FU Orionis observed outburst ([Clarke et al.l 
Il99d ; iBell fc Linl 1 19941 ; lLodato fc Clarkd I2004T ). ass» 
0.001 — 0.003. In fully ionized discs in dwarf novae, the best 
observed evidence suggests a typical ass ~ 0.1 — 0.4, whilst 
the relevant n umerical MHD sim ulations evaluate ass ~ 10 
times smaller. iKing et al.l (|2007f ) attribute this discrepancy 
to incorrect magnetic and boundary layer shortcomings in 
the computations. Nevert heless, this i s not the co nclusion 
of the full story because ILanzafamd (|2008l . I2009T I showed 
that a well bound viscous accretion disc structure mod- 
elling strongly depends on several conditions: the kinematic 
of the mass transfer, 7, ass and so on. For isothermal or for 
a quasi-isothermal thermodynamics, a disc is structurally 
bound even for ass = because numerical dissipation only 
is enough to produce a disc in s h ocks events, a resul t also 
discussed bv lSawada et al. |l987l); Spruit et all (|l987h and, 
more recently, by ILanzafamd ( 2010a ) in its physical sense. 



4 v AND C COEFFICIENTS AND 

MOLECULAR CHARACTERISTICS 

Different formulations of macroscopic physical dissipation 
rarely converge with each other and often include an arbi- 
trary parameter, to be evaluated. In addition, any corre- 
lation to microscopic (molecular, atomic, nuclear) physical 
properties is absent. 

To this purpose, we propose the following evaluation of 
the physical kinematic viscosity coefficient v, to be used to 
determine the viscous stress tensor r in the Navier-Stokes 
equations. 

Microscopic molecules, atoms, nuclei, have known elas- 
tic scattering impact cross section k, useful to compute v. 



Without any consideration to the existence or to the involve- 
ment of internal energy levels for an ideal flow, physical dis- 
sipation transfers macroscopic ordered kinetic energy flows 
into heat, that is in microscopic chaotic kinetic energy flows. 
This means that elastic scattering collisional cross sections 
for an ideal gas are those right for a reformulation of v. For a 
gas mixture, the mean value of the elastic scattering impact 
cross section 



(23) 



should be considered, where Xi — Ti»/^. rii is the rel- 
ative numerical abundance of the chemical species i. 

However, we need to take into account the total number 
of microscopic molecules-atoms within the macroscopic cross 
section determined by the mixing length as in eq. (18) which 
arithmetically gives as a result a mixing length close to the 
smallest scale length in the summation. For this reason, we 
prefer to compute I as 



I = min(ii, h, h, — , In), 



(24) 



\v\/V 



the various U referring to p/|Vp|, p/\Vp\, 
|u|/|V x v\, and so on. 

Hence, the total number of barions contained within the 
3D mixing mass pi 3 is pl s /jlmH, where mn is the proton 
mass and Jl is the mean molecular weight. 

Since we need the molecular collisional cross section, 
in 3D it is necessary to compute (pi 3 /Jltuh) 2 ^ 3 , that is the 
number of molecules within the volume Z 3 , powered to 2/3. 
This number has to be multiplied with k to get a statisti- 
cally effective collisional surface composed of a multitude of 
microscopic cross sections: 



f 

pm H 



2/3 



(25) 



To calculate u, we need to divide this arithmetic term 
by a length A. This length decreases whenever the local nu- 
merical density increases. Therefore A ~ n -1 ' 3 . 

As far as the velocity contribution of v is concerned, 
we exclude not only Vfi ow , but also the thermal e 1//2 , being 
this last strictly linked to the chaotic microscopic kinemat- 
ics. Kinematic velocities of extraneous bodies - as in the 
original Von Karman formulation - moving in the fluid are 
not considered. This exclusion is a consequence of the fact 
that the kinematic viscosity coefficient, mutually correlated 
to the thermal conductivity and the diffusivity coefficients, 
are all expressions of the intrinsic physical property of the 
fluid. Hence, to conclude, our 



P i 3 



pm H 



2/3 



— 1/3 

nn c s 



Pi 2 ^ 



pm H 



(26) 



being n = p/jlmn- 

Its 2D counterpart, according to the same algebraic log- 
ical steps is: 



pm H 



1/2 



1/2 

KTl C s = 



= 7 

(im H 



(27) 



being E the 2D mass density. 

In these expressions, both the molecular/atomic Jlmn 
and k are included, as well as other macroscopic physical 
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quantities like E or p. This means that, the spatial com- 
ponent of v could strictly be not shorter than I and longer 
than h, according to other physical variables, especially E 
or p. This formulation is free of any arbitrary parameter. 
Moreover, in (26) and (27), the quantities (p /jEtTn,Er)iH and 
(S//Img)K are pure numbers both < 1, otherwise the den- 
sity of the fluid is comparable or greater than the atomic 
density. It is important to note that eqs. 26 and 27 are not 
exactly equivalent. In eq. 26 v oc I 2 in its spatial dependence 
(as Prandtl's eqs. 16, 17), while in eq. 27 v oc I in the same 
spatial contribution (as Shakura and Sunyaev's eq. 22). Eqs. 
26 and 27 should be considered strictly correlated either to 
a 3D or to a 2D modelling, respectively. Indeed, densities p 
and E are formally correlated by the equivalence of their nu- 
merical densities as (E/7im_f/) 3 = {p /JIitih) 2 ■ So that E = pi 
would be a very specific case. 

Molecular or atomic collisional cross sections are as- 
sumed, for the sake of simplicity, circular, without any dis- 
tortion. Molecular dimensions are determined by the so 
called Van der Waals mean radius, defining the limit where 
the microscopic molecular force potential becomes attrac- 
tive, deviating from that of a non interacting free particle 
relative to an ideal gas. In the case of fully ionized gas, the 
repulsive Coulomb elastic scattering among head on collid- 
ing ions determines the shortest classical impact parameter 
r p as r p ~ (ATV€ )~ 1 2ZiZ2e 2 /SKbT, related to k as k = irr 2 ,. 
Z\e and Zz£ are the two effective electric charges of the two 
colliding ions, Kb is the Boltzman constant, T is the tem- 
perature and e Q is the dielectric constant of the vacuum. 

Notice that the decrease of the mixing length I in a more 
effective characteristic scale length as a result of the presence 
of a multitude of small scale constituents, as shown in eqs. 
(26, 27), does not alter the meaning of the kinematic viscos- 
ity coefficient role in the tensorial expression of the stress 
viscous tensor (eqs. 12, 13) in the Navier-Stokes equations 
describing fluid flows. Therefore, although we introduce lo- 
cal physical properties of the fluid, and although turbulence 
is not a feature of fluids but of fluid flows, the reformulation 
of the scalar v within T a ,p stays meaningful. 



4.1 Kinematic viscosity and thermal conductivity 
coefficients for dilute gases 

The kinematic viscosity coefficient v and the thermal con- 
ductivity coefficient c are dimensionally identical. Both are 
characterized by a scale length multiplied by a scale veloc- 
ity. Hence, both are transport coefficients. A relevant dif- 
ference between viscosity and thermal conductivity is that 
while viscosity is activated only whenever a relative mo- 
tion occurs among contiguous flow elements, thermal con- 
ductivity is an energy transport mechanism always existing 
whenever a temperature gradient occurs even in steady state 
conditions. However, the two coefficients are always related 
to each other because both explains the tendency of the 
thermodynamic system toward a homogeneity and isotropy 
status, smoothing out kinematic (y) a nd therma l (c) local 
spatial discrepancies. The ratio c/v |R]|eif (1 965)) is: 



cy 



(28) 



c v = 3i?T/2 = 3K B T/m H - Experimentally {c/v){ep/c v ) « 
1.3 — 2.5, instead of 1. The discrepancy is largely explained 
because of the fact that theoretically c is evaluated consider- 
ing a uniform molecular velocity distribution instead of local 
molecular kinematic differences related to the presence of a 
temperature spatial gradient. As a consequence, considering 
v = £, for an ideal gas, 

K B T 



tji ep epran 
so that, 

c ^ 2 ( I V p^c s «L = 3K B T ( J-) 2 p~R— 
\p J m H e \pm H J e 

in 3D, and 



i y a jl Cs 2k = 3KbT ( j_ y E£«— 

p J m H e \ pm H I e 



(29) 



(30) 



(31) 



where cv is the molar specific heat of the gas at 
constant volume which, for an ideal (monoatomic) gas is 



in 2D. Since coefficients v and c are comparable for 
dilute gases, then the viscous and the thermal conductivity 
time scales are also comparable. This means that the inertia 
of matter tends toward a uniform kinematic and thermal 
configuration with the same time scales. 

In the viscous computational fluid dynamics, currently 
v and c are not only arbitrarily parametrized, but also not 
correlated from each other, as it correctly should be. In spite 
of the fact that results in isothermal or in nearly isothermal 
conditions could still be significant, the lack of any correla- 
tion between v and c is free of any physical meaning. 

In the rest of the paper, we mainly pay attention to the 
physical viscosity. However, since now onwards, any con- 
clusion referring to the role of the physical viscosity will 
also refer to the role of the thermal conductivity in a close 
cause-effect correlation, where any any significant local spa- 
tial derivative involving a relative motion among contiguous 
flow parts will involve a a braking and a viscous heating; any 
local heating will involve larger pressure and temperature 
gradients and consequently a transport of energy and mass 
with comparable time scales; a heat transfer will involve a 
decrease of temperature and pressure spatial gradients. 



5 VISCOSITY TESTS 

Flow transport and damping of turbulence are the only two 
fields where the role of physical dissipation in the Navier- 
Stokes equations is better emphasized. Therefore, in this 
section, appropriate tests on the physical viscosity efficiency 
will be done. In particular, regard a 2D SPH shockless mod- 
elling on the radial flow viscous transport in an annulus ring 
and a 2D SPH modelling of Burger's turbulence. 7 = 5/3 will 
be assumed throughout the tests. This low compressibility 
regime is the less advantageous condition for viscous forces 
against pr essure forces, as well as against numerical-artificial 
damp ing l|Molteni et al.lll99ll ; iMurravl [l99l lOkazaki etafl 
2002), explained by a modest contribution of the bulk com- 
ponent V-w in eqs. (12, 13), not compensated by the high c s , 
on v. In such tests, results on v — f , also including a thermal 
conductivity cot v will also be compared with those relative 
to 



c s h 



(32) 
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v = c s l (33) 

which are the two simplest analytical expressions for 
v, where the characteristic length could be either at its 
minimum geometric value h, or coincident with the mix- 
ing length I. Despite their full mathematical meaning, for- 
mulations (32) and (33), like (eqs. 16, 17, 22), lack of any 
physical sense without any correlation to what matter the 
flow is made of. 

5.1 2D radial viscous transport in a shockless 
isothermal annulus ring 

Theory on 2D shockless radial transport in a Keplerian an - 
nulus ring in a gravitational potential well (jPringld Il98ll ) 
predicts that, from the Green function, the solution of the 
initial Keplerian mass distribution at time t = for E is: 

E(r, t = 0) = mS(r - r )/2nr (34) 

for a ring having mass m and an initial radius r . The 
solution, at time t, in terms of dimensionless radius x — r/r a 
and viscous time 9 = Vlvtr^ 2 is 

E(x,t) = (m/7rr~ 2 )6|- 1 a;- 1/4 exp[-(l+a; 2 )/6»]J 1/ 4(23;/6»),(35) 

where I\ia is the modified Bessel function. The action 
of viscosity is to spread out the entire annulus ring toward 
a disc structure transporting most of the low angular mo- 
mentum mass toward the centre of the potential well and 
transporting a smaller fraction of high angular momentum 
mass toward the empty external space. 

For practical computational purposes, since it is impos- 
sible to reproduce a delta Dirac function at time t = 0, it 
is necessary to start numerical calculations from an initial 
mass distribution relative to a small 9 value. In our example, 
we choose 9 = 0.017, a value comparable with that used by 
ISpeith fc Klevl l|2003h for a radial viscous transport similar 
test. 

Fig. 1 shows a comparison among several radial distri- 
butions of the 2D mass density E(r, 9) for various 9. Because 
of the universality of graphical representation by using 8 in- 
stead of t, computational radial profiles at an assigned 9 
value have to correspond with each other, fitting the ana- 
lytical solution (eq. 35). In particular, fig. 1 shows in the 
same picture E(r, 9) profiles both for two Shakura and Sun- 
yaev kinematic viscosity v — assCsH, and for the v = c s h, 
where h — 0.09, and for the newest v = f , where our for- 
mulation for v is used. This clearly shows that the viscous 
v = £ model correctly replicates the right radial shockless 
transport mechanism, without any local distortion. In this 
example, we used throughout the models, M a = 2 • 10 33 g, 
Ro = 10 11 cm, as typical astrophysical values. We considered 
an initial density E of the order of 10~ 10 g cm~ 2 , an initial 
sound velocity c s — 5 • 10~ 2 « o , where v = 2n(GM /R ) 1 ^ 2 
cm s" 1 is the normalization value for the velocity, and a gas 
composed of pure molecular hydrogen. The thermal energy 
per unit mass e is kept constant, being de/dt — throughout 
the entire simulation. Variations of these parameters do not 
produce any difference in the radial density distribution for 
each 9, being the viscous time 9 an absolute reference time 
for E. The v — c s h profile is shown as a further test on the 
radial profile, where eq. (32) has been considered in the v 
calculation. Instead, models adopting eq. (33), or even the 
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Figure 1. Examples of four calculated mass radial distributions 
as a function of the viscous time 8. Two a Shakura and Sunyacv 
radial distributions are shown, as well as a v = c s h, with h = 0.09, 
compared with our v = § radial profile. 

Prandtl formulation (15), do not yield any realistic radial 
profile of mass distribution because of the too large viscos- 
ity, due to the very negligible spatial derivatives in the I 
calculation, causing the rapid accretion of the entire disc. 

In spite of the low spatial resolution adopted for prac- 
tical purposes, results on radial shockless viscous transport 
here reported clearly show that, without any other physical 
alternative, the Shakura and Sunyaev formulation looks like 
an appropriate expression for v for 2D shockless disc struc- 
tures, since it uses a characteristic length H much shorter 
than I along a 3rd dimension that does not exist in 2D. How- 
ever, physical formulations as in eq. (26) or (27), where the 
effective scale length (Ei/7imH)^ is much shorter than / in 
the case of diffuse matter, could be a valid physical alterna- 
tive, free of any arbitrary parameter. 

5.2 Damping of 2D Burger's turbulent flow 

Statistical studies of turbulence normally involve hypothe- 
ses about homogeneity and is otropy on the distri bution and 
kinematics of the spatial flow ( Kol mogorovll 1941al lbl) . 2 D tur- 
bulence is relevant to understa nd large scale flows l|Frishl 
1 19951 ; iKellav fc Goldburd |2002| ). 2D turbulence schemati- 
cally discusses either a "forced steady state turbulence", or 
a " decaying turbulence" , if an explicit forcing term is added 
in the momentum equation or not. Eddies of different sizes 
showing density and potential fluctuations in the flow nor- 
mally characterize turbulence. 

2D phenomenology is somewhat more complex than 3D 
pheno menology, even though computationally more conve- 
nient (|Tabelinel [2002). being not derivable from simple di- 
mensional arguments. 3D turbulence involves scales smaller 
than the trigger one and it is supposed to be hosted to a 
direct enstrophy cascade from large to small scales (from 
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Figure 2. (X, Y) plots of density map at various times T for the 
physically non viscous model with h = 0.05. 64 greytone are used. 
»x, vy tomograms are also reported, showing velocity fluctuation 
both during the initial turbulent phase and during the following 
damping subsequent phase. 



small to large wavenumbers), where the mean kinetic en- 
ergy is transferred and mean enstrophy is conserved. In- 
stead, 2D turbulence involves larger scales than the trigger 
one and it is supposed to be hosted to an inverse enstro- 
phy cascade from small to large scales (from large to small 
wavenumbers). An inverse cascade of energy and a contem- 
porary d irect cascade of en strophy in 2D is called a "dual" 
cascade IjManz et alj 12003 ). Even though a physical viscos- 
ity is usually considered, as responsible of flow damping of 
the Navier-Stokes equations, the phy sical counterpart of ar - 
tificial viscosity has been discussed in iLeVeaudll 19921 . 120021 ); 
iFletcherl l|l99lf ); Iffirscbl i|l997l ). 



Figure 3. (X, Y) plots as those of Fig. 2 for the physically viscous 
model v = c s h. 



Turbulence cascade theory substantially predicts that 
the kinetic energy is contained in turbulent eddies, while 
the cascade process for enstrophy is somewhat ambivalent in 
2D, according to boundary and initial conditions. The larger 
the eddy, the larger its kinetic energy content, according to a 
power scaling law of eddy dimension, respecting the energy 
conservation law (kinetic + thermal) within the whole sys- 
tem. In 2D, the so called "dual cascade" process determines 
the formation of larger turbulent eddies up to the limit of 
the entire spatial domain. 

The numerical experiment here carried out for 2D 
Burger's turbulence studies the temporal evolution of the 
damping of a chaotic gas inside a L = 5 size 2D squared 
box, whose density and whose kinetic velocity are initially 
locally random, as an example of decaying turbulence, vx 
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Figure 4. (X, y) plots as those of Fig. 2 for the physically viscous 



model v ~ « 2 ||2i I 
I m 



Figure 5. (X, y) plots as those of Fig. 2 for the physically viscous 
model v = c s i. 



and vy values are within the range —5 • 10 -4 and 5 • 1(T 4 at 
the beginning while h = 0.05. Hence vx ~ and vy ~ sta- 
tistically at time T = 0. Throughout the models the initial 
adimensional thermal energy per unit mass e = 1.86 • 10~ 7 . 

The initial settings consists of an uniform thermal en- 
ergy together with a statistical macroscopic homogeneous 
and isotropic spatial distributions on density and on veloc- 
ity components vx, vy- Being the macroscopic scale length 
much larger than the spatial resolving power, clear fluctua- 
tion exists in the homogeneity on the scale resolution length. 
This is enough to ignite the initial turbulent kinematics, 
with larger velocity, thermal energy and density fluctuations 
as larger are pressure forces from the beginning. Damping 
viscous effects always work since the beginning. Hence, as 
a decaying test, after the initial turbulent condition, where 



turbulence dominates, the final configuration is character- 
ized by a fluid statistically at rest, where fluctuations in the 
velocity field, in the thermal energy and in the density are 
reduced after some time toward a more general uniformity. 

Figg. 2 to 6 show XY distribution of density at various 
times T for various physical kinematic turbulent viscosities 
v, as well as parallel velocity tomograms showing vx and vy 
at the same time, from time T = to T ~ 100. 64 grey tones 
for the density maps are shown between the minimum and 
the maximum E values for each plot. That is, each plot has 
its own minimum and maximum value for E. The normaliza- 
tion value for the density is E = 10~ 10 g cm~ 2 and densities 
are initially computed multiplying E times a random num- 
ber between and 1. We distinguish two physical regimes. In 
the first one, from the beginning to T « 25, 2D turbulence 
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Figure 6. (X, y) plots as those of Fig. 2 for the physically viscous 
model v = £ for a H2 gas. S Q = 10~ 10 g cm~ 2 . 



Figure 7. (X, Y) plots as those of Fig. 2 for the physically viscous 
model v = £ for a O2 gas. So = 1.6 ■ 10 -9 g cm -2 . 



dominates because small eddies grow up. Linear size of ed- 
dies, initially of the order of h, grows up to several h values. 
The relevance of the artificial dissipation - v m (0.1 — 0.4)c s /i 
- is evident in the physically non viscous case, because of the 
low spatial resolution adopted. From T ~ 25 onwards, clear 
differences both on the density spatial distribution and on 
the velocity component fluctuation appear, distinguishing 
the efficacy of the various adopted physical viscosity coef- 
ficients, since the physical viscosity better works increasing 
c s , as well as enlarging the eddies if v oc I or oc I 2 . Plots 
clearly show that v — c s l is the most dissipative model in 
so far as E = 10 -10 g cm~ 2 . In this example, the physical 
damping relative to v — £ is comparable with that relative 
to the Prandtl's v ~ l 2 \dvdx\. 

A reduction of the intrinsic damping by decreasing the 



spatial resolution length is possible, of course, but not prac- 
tical. The reduction of the time step, together with a much 
larger computer memory needed involve very long compu- 
tational time, even for 2D simulations. We adopted a low 
spatial resolution, since h/L — 1CP 2 . This involves that re- 
sults inherent to all models are quite viscous. However, these 
results, even though conditioned by the intrinsic numerical 
damping effect, clearly show how £ correctly works, free of 
any dependence on h, or on any arbitrary parameter. 

To show the role of other chemical species, through their 
mean molecular weight, we also show in Fig. 7 results regard- 
ing the damping of 2D turbulence for a gas of molecular O2 
for E = 1.6- 1CP 9 g cm~ 2 at time T = 0. These results have 
to be compared with those for a gas of molecular H2 with the 
same initial numerical density (Fig. 6), whose E = 10 -10 
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Figure 8. (X, y) plots as those of Fig. 2 for the physically viscous 
model v = £ for a i/2 gas. So = 10 _s g cm -2 . 



g cm - at T = 0. Results for O2 are influenced by the gas 
molecular cross section 75, in so far as the numerical density 
p/jlmH is relevant in the £ calculation to get a v comparable 
or higher than any artificial or any numerical dissipation. In 
this example, being the numerical densities for the molecu- 
lar O2 at the beginning the same as those of Fig. 6 for H2, 
v — £ is so effective that viscous dissipation quickly damp- 
ens both any initial velocity fluctuation and flow progression 
toward any spatial homogeneity. 

An even stronger viscous flow damping is shown in Fig. 
8 for a denser gas of molecular H2, whose E = 10~ 8 g cm~ 2 
at T = 0. Viscous molecular damping, where v = £, appears 
more effective for higher density flow modelling. As an ex- 
ample, in Fig. 8 viscous dissipation is so relevant that any 
local flow motion is strongly suppressed since the beginning, 



also because of a more relevant heat thermal conductivity 
suppressing local thermal differences caused by the viscous 
heating. At the same time, the viscous heating is so strong 
that the consequent increase of pressure forces time by time 
activates again the velocity fluctuations in the flow after 
T ~ 50. Hence, a different evolution of velocity fluctuation 
is developed, depending on the initial conditions for v, e, v, 
E. The study of the 2D viscous damping of the 2D kinemat- 
ics of a so dense fluid does not need to concern us. 



6 3D ACCRETION DISC IN A CLOSE BINARY 

We compare results for 3D stationary disc structures both 
in high compressibility (7 = 1.01) and in low compressibil- 
ity (7 — 5/3) with the aim of getting a physically well- 
bound accretion disc around the compact primary star in 
a LMCB. However, these comparisons are thought to show 
which physically viscous disc structure, in the Shakura and 
Sunyaev formulation (eq. 22), in the Prandtl formulation 
(eqs. 16, 17), and in the more simple formulations (eqs. 32 
and 33), better compare with that relative to v = £. The 
adopted spatial resolution length is h = 5 • 10~ 3 through- 
out the simulations, being the mutual separation of the two 
stars normalized to 1. 

The characteristics of the binary system are determined 
by the masses of the primary compact white dwarf and of 
its companion star and their separation. We chose to model 
a system in which the mass Mi of the primary compact star 
and the mass M2 of the secondary subgiant star are both 
equal to IMq and their mutual separation is di2 = 10 6 Km. 
The injection gas velocity at LI is fixed to Vi„j ~ 130 Km 
s _1 while the injection gas temperature at LI is fixed to 
To = 10 4 K, taking into account, as a first approxima- 
tion, the radiative heating of the secondary surface due 
to radiation coming from the di sc. Superso n ic kinemati c 
conditions at LI are di scuss ed in iLanzafamel (|200Sl 120091 1; 
lLanzafame et all (|2000l . l200ll ). especially when active phases 
of CB's are considered. The reference frame is centred on the 
primary compact star, whose rotational period, normalized 
to 27r, coincides with the orbital period of the binary system, 
being the velocities normalized to v = [G(M\+M2/d\2)] 1 ^ 2 ■ 
Results of this paper are to be considered a useful test to 
check whether disc structures (viscous and non) show the 
expected behaviour. 

We simulated the physical conditions at the inner and 
at the outer edges as follows: 

a) inner edge: 

the free inflow condition is realized by zeroing gas flow inside 
a sphere of radius 10~ 2 , centred on the primary compact 
star. Although disc structure and dynamics are altered near 
the inner edge, these alterations are relatively small because 
they are balanced by a high gas concentration close to the 
inner edge in supersonic injection models. 

b) outer edge: 

gas flow from LI towards the interior of the primary Roche 
Lobe is simulated by constant gas pressure, density, and 
thermal energy per unit mass, as well as a constant veloc- 
ity in a small conic region having LI as a vertex and an 
aperture of ~ 57° . The radial length of this small volume is 
~ Wh. The initial injection particle velocity is radial with 
respect to LI. Local density at the inner Lagrangian point 
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Figure 9. (X, Y) plots for 64 greytones of density map for 3D physically viscous accretion disc in the v = cess c sH Shakura and Sunyacv 
formulation for 7 = 1.01. Time T and ass are a l so shown. 



LI: pL\ = 10 -14 g cm -3 . This order of magnitude is ex- 
plained as follows: from the flux conservation, pvS = const. 
Hence, piviSi — P2V2S2 from the two Roche lobe sides at 
LI. If the two stars have comparable masses, piV\ ~ p2V2- If, 
at LI, from the secondary subgiant star side, p2 ~ 1CP 11 g 
cm -3 and vi ~ 10 _1 Km s -1 as typical photospheric values 



for a subgiant, then pi 

-1 



10" 



g cm , being Wi w 10 Km 



Supersonic ma ss transfer conditions fro m LI were pre- 
viously adopted in Lanzafamc (2008, 2009), where disc in- 
stabilities, responsible for disc active phases of CB are dis- 
cussed in the light of local thermodynamics. Whenever a 
relevant discrepancy exists in the mass density across the 
inner Lagrangian point LI between the two stellar Roche 
lobes, a supersonic mass transfer occurs as a consequence of 
the moment um flux conservation . The same result can also 
be obtained (|Lubow fc Shull 19751 ) by considering either the 
restricted problem of three bodies in terms of the Jacobi 
constant or the Bernoulli's theorem. Although the stellar 
gas is surely neutral for T = 10 4 K at LI, we consider a 



chemical composition of pure ionized hydrogen for the sake 
of simplicity for a simple calculation of k in the disc bulk. 

Figg. 9 and 10 show five XY plots of the physically vis- 
cous discs for 7 = 1.01 and for 7 = 5/3, respectively, in 
the Shakura and Sunyaev prescription, where ass ranges 
from 0.1 to 0.5 in steps of 0.1. Each density map has its 
own minimum and maximum densities, scaled in 64 grey- 
tones. These ass values are in accordance with the typical 
and with the maximum compa tible with both astrophysical 
observa tions jKing et all 2007) and with numerical exper- 
iments (|Abolmasov fc Shakura! [2009) for these kind of as- 
trophysical objects. As it is evident, this formulation always 
works quite well, complicating the determination of the right 
value for the ass arbitrary parameter up to the point that 
viscous disc structures are practically indistinguishable from 
each other for 7 = 5/3. 

Figg. 11 and 12 also compare density maps for 7 = 1.01 
and for 7 = 5/3, respectively, both for the physically inviscid 
and for the other four physically viscous formulations for v. 

In the non viscous model, for 7 = 1.01, Fig. 11 shows 
a compact disc structure, where a high density central 
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torus is formed in the disc bulk. In this case, the intrinsic- 
artificial viscosity alone favours the radial mass viscous 
transport, as well as the conversion of mechanical energy 
into heat. A glimmer of spiral shape structures at the 
disc outer edge exists, although not well developed because 
of the limited radial extension of the disc together with 
its substantial circular symmetry. The tight link between 
the elliptical extension of eccentric disc and its asymme- 
try to the presence o f spirals coming from it s out e r edg e 
has been d i scusse d by iBisikalo et ail (|l998al lbl. Il999l . l2000h ; 
lLanzafamei (|200jK especially in relation to the "tidal trun- 
cation radius" dPapaloizou fc Pringlel Il977l ; IZang fe Chenl 
Il992l ; llchikawa fc Osakilll992l . ll994 V 

Notice in Fig. 11 for 7 = 1.01 that, for a gas of pure 
ionized hydrogen, the viscous disc structure whose v = £ 
tightly compares with that whose v = assc s H, where 
ess = 0.3 — 0.5 (in Fig. 9). In the same figure 11, physi- 
cally viscous discs, whose v ~ l 2 \dv/dx\, v = c s h and v — £, 
are also elliptically extended, with a clear evidence of spiral 
structures coming from the eccentric disc outer edge from its 



more lengthened zone and from the disc side where the flow 
stream coming from LI collides to the disc outer edge. For 
v ~ l 2 \dv/dx\, the existence both of a disc thickness and of 
collisional spatial derivatives limits I, and consequently v, to 
local values not so large as those relative to the 2D shockless 
radial viscous transport for an annulus ring in §5.1. 

In 3D, despite the better limitation for Z, the tendency 
toward a larger dissipation, as previously seen in §5.2 for 
the damping of Burger's turbulence for 2D structures (or 
for flattened 3D structures), gives the v = c s l disc struc- 
ture a quite viscous behaviour. In the highly compressible 
and highly viscous case, not only the disc radial extension 
enlarges, but at the same time the entire structure bet- 
ter shows a retrograde precession in a waving of its outer 
edge, according to the reference frame of Fig. 13. This ret- 
rograde precession is originated by the gravitational tidal 
forces and supported by a significant Coriolis acceleration 
caused by the strong shockless radial viscous transport in 
a spiral-shape kinematics. Notice how in the momentum 
equation, the Coriolis acceleration is the only term depend- 
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Figure 11. (X,Y) plots for 64 greytones of density map for 3D physically viscous accretion disc for 7 = 1.01. Time T as well as the 
used v formulation are also shown. 



ing on the radial component of v. The role of the Corio- 
lis acceleration in the exchan ge of fluid blobs in disc struc- 
tures is well described in the lFrank et alj |2III)2) textbook. 
The high viscosity and the high compressibility work to- 
gether producing a gas structure comparable with that of a 
soft rubbish elastic membrane. Being the fluid highly de- 
formable, but at the same time structurally bound, it is 
highly characterized by uniform collective motions and by 
an efficient conductive thermal transport, where any per- 
turbation, even though strongly attenuated, propagates and 
involves as many fluid parts as possible. As a consequence, 
disc outer edge distortions are much more stressed. The disc 
distortion and the consequent lengthening of its outer edge 
better enhances the effect of tidal forces on the entire disc 
fluid dynamics. Results showing a retrograde clockwise pre- 
cession in the viscous disc outer edge in a longitudinal os- 
cillation dis t orting th e entire disc profi l e were obtained by 
IWhitehurstl |l 988al lbh ; iKlev et al. et"ai] (|200ST l, also in high 
artificial viscosity conditions, for LMCBs. Hence, a high vis- 
cosity together with a high compressibility gas are the es- 
sential conditions to get accretion discs where disc outer 



edge distortion and its tidal precession are evident in op- 
position to the standard disc models where such dynamics 
is absent. For v — c 3 l disc model, what is anomalous is 
the short period of waving of its distorted outer edge in 
a LMCB. It looks like comparable with the orbital period. 
Whenever in a LMCB a retrograde precession of the disc 
outer edge exi sts, it is at least fO times longer than th e 
orbital period l|Whitehurstlll988al lb1; [Klev et al. et alj|2008h . 
especially if tidal forces activate both the disc warping and 



its outer edge precession jKatdl 19731; Iwf icrs fc Prinridfl999l ; 
iBate et al.l l200d ; lOgilvie fc DubuslhoOOT ). This anomaly is 



however simply explained considering that for v = c s l the 
viscous time scale ~ r 2 jv is much shorter with respect to 
the viscous time scales of other less viscous disc modelling. 

Disc phenomenology is different in the low compressibil- 
ity 7 = 5/3 modelling, as shown by the comparison between 
Figg. 10 and 12. In the non viscous case, pressure forces as 
well as mass outflow from the disc outer edge are so rele- 
vant, that the entire disc structure is sparse and not well 
bound within the primary's gravitational po tential well of 
a LMCB, a consolidated result since works of iMolteni et alj 
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Figure 12. (X, y) plots for 64 greytones of density map for 3D physically viscous accretion disc for 7 = 5/3. Time T as well as the used 
v formulation are also shown. 



|l99ll ); lLanzafame et al.1 l| 19921 ). The modest physical vis- 
cosity from the Prandtl's and from v = c s h does not allow 
any binding of the flow into the primary's gravitational po- 
tential well. Only a tiny gas concentration toward the inner 
regions of the potential well is visible. Instead, the two disc 
models for v = £ and for v — c s l look like working quite 
well, being also in a good comparison with the low com- 
pressibility viscous discs in the Shakura and Sunyaev pre- 
scription. In these two highly viscous collisional disc mod- 
els, the larger physical viscosity damping, also evidenced in 
the 2D viscous damping of Burger's turbulent test, as well 
as an efficient thermal conduction smoothing out thermal 
and pressure gradients, determines an effective disc binding 
within the gravitational potential well. Moreover, the rele- 
vant high pressure forces for 7 = 5/3 prevent any amplified 
disc oscillation of its outer edge as previously evidenced in 
Fig. 11 for 7 = 1.01 and v — c s l. 



7 CONCLUDING REMARKS 

In this paper we formulate a physical general expression for 
the kinematic viscosity coefficient v, able to determine a cor- 
rect physical viscous flows in the nonlinear Navier-Stokes 
fluid dynamics. In this formulation, we search for an ex- 
pression, free of any arbitrary numerical parameter, pay- 
ing attention to the correlation between the microscopic 
molecular-atomic cross section to macroscopic characteris- 
tic lengths (mixing length) related to the linear eddy di- 
mension. At the same time, a reformulation for v in more 
strictly physical terms also involves a reformulation for the 
thermal conductivity coefficient c, being v and c simply cor- 
related for a dilute gas. 

Current adopted kinematic viscosity coefficients are 
normally written as a pure mathematical coefficients, with- 
out any correlations to the molecular cross sections. Some- 
times, as for example in the Shakura and Sunyaev formula- 
tion, limited only to disc geometric structures of the flow, 
an unknown arbitrary parameter appears. 

Current v formulations are not used whatever is the 
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Figure 13. (X, Y) plots for 64 greytones of density map for 5 instants of the 3D the physically viscous accretion disc for 7 = 1.01, 
v = c s l. Time T is also shown. 



viscous physical problem. Either v is specific for shockless 
viscous transport phenomena, or it is specific for chaotic 
turbulent flows. Of course in both cases, the conversion of 
mechanical energy into heat, together with the braking kine- 
matics, yield laminar flows after some time without any ex- 
ternal force acting on the flow. 

fn this paper we also pay attention to the role of the 
intrinsic flow damping. It could be either explicit, as an 
artificial viscosity term, or it could be numerical, or both, 
especially for finite difference codes, because of second or- 
der terms by the numerical conversion of spatial derivatives 
coming from the Taylor series expansion. This non physical 
damping is necessary to handle collisional shocks, prevent- 
ing any distortion in the shock fronts. In both situations, the 
adopted spatial resolution h has a direct or indirect role. 

For practical reasons, we worked in a low spatial resolu- 
tion, adopting h/L — W~ — 10~ throughout. This means 
that our results are affected by a relevant intrinsic damp- 
ing. However, in spite of this disadvantage, we clearly stated 
that: 



• the v — £ correctly determines a physical radial mass 
viscous transport in an annulus ring, as it is in the case of 
the Shakura and Sunyaev formulation. If the characteristic 
mixing length, determined by expressions like (17) or (23), is 
comparable with L, other formulations for v oc I, or v cx I 2 
(eqs. 16, 17, 33) looks like too viscous, so much that the 
entire disc is quickly accreted. 

• the role of the numerical density, as well as both of the 
molecular-atomic cross section and of eddy linear dimension 
(I) on v — £ are evident in so far as the physical damping 
is at least comparable with that of the intrinsic dissipation. 
Eqs. (26 and 27) include all physical situations. The main 
contribution on v — £ either comes out from a long I, as in 
the case of the mass radial viscous transport in the annulus 
and in disc structures, or it comes out from a large numerical 
density, as in the case of random chaotic motions. Neverthe- 
less, specifying the chemical composition of the fluid, what 
is relevant is not only I, as eqs. 26 and 27 show. 

• Expressions where v oc h are arbitrary, of course, de- 
spite low resolution results ofd simulations are physically 
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meaningful, because there is not any physical correlation 
between the kinematic viscosity and the eddy linear size. 

• In the 3D accretion disc simulations, results on steady 
disc structure for v = £ are in a good accordance to those rel- 
ative to the Shakura and Sunyaev formulation without any 
ambiguity on the arbitrary ass parameter. The comparison 
in both high and in low compressibility are very success- 
ful. In particular, for 7 = 1.01, v = £ looks like working as 
v = assc s H with ass = 0.3 — 0.5. Instead, other expressions 
for v arc unsuccessful either in one case or in the other case 
in so far as the Shakura and Sunyaev viscous prescription is 
correct. On the contrary, the highly viscous disc modelling 
could deeply pay attention to the relevance of the v = c s l 
kinematic viscosity coefficient in the solutions of the Navier- 
Stokes equations. It is true that high viscosity disc models 
could also be built up tuning either p and/or k for v = £. 
However, high viscosity physics could be obtained in high 
density conditions (up to terrestrial conditions, see eqs. 26, 
27) that very rarely happens in diffuse matter astrophysical 
environments. 

For astrophysical objects, showing flat accretion disc 
structures, the v — assc 3 H Shakura and Sunyaev formula- 
tion is successfully adopted. In this case, H m TdiscCs/v^, « 
10~ 2 r disc for a correct application in a strict Keplerian tan- 
gential kinematics, ass ranges within ~ 0.001 — 0.4 accord- 
ing to the astrophysical object considered, not without any 
ambiguity. The only physical local component in this formu- 
lation is the sound velocity, that is assumed as a function of 
the mass accretor and of the radial distance only, in shock- 
less conditions. It does not exist any consideration on why 
the Prandtl's v ~ l 2 dv/d x is not taken into account. In 
our formulation, every component of v — £ depends on the 
local conditions, with or without shock events. This, with- 
out any doubt complicates an exact calculation on v since, 
apart the chemical composition, the entire (£/7i?Tiij)Dt or 
(p/jZrriH)l k terms also depend on the local densities and 
on I, which could be very different from the local disc thick- 
ness H. This without any consideration on k that could be 
very different from molecular-atomic scales to high temper- 
ature nuclear scales of cosmological interest. As an example, 
for AGN, the term H/jimH w 0.1 — 1 g _1 cm 2 for protons in 
so far as T « 10 6 — 10 7 K. This means that if we arbitrarily 
impose I = H, either S or pi « 10~ 2 g cm~ 2 . Being the first 
2D E value too high, instead the second 3D pi value looks 
like much more plausible. For star forming objects like FU 
Orionis, the ass ~ 10 -3 in the v = assc s H Shakura and 
Sunyaev formulation. If we consider a gas of pure neutral 
atomic hydrogen, the ratio ~K/~pniH ~ 10 4 g _1 cm~ 2 . Hence, 
S or pi ~ 10~ 12 g cm~ 2 , that means a quite high value for 
S or a low, but more realistic value for pi. The entire eval- 
uation on v = £ so far includes molecular-atomic-nuclear 
characteristics only on the basis of a collisional molecular 
gas dynamics without any consideration on the role of elec- 
tric and/or magnetic effects, as well as on the presence of 
dust in the diffuse matter. On the contrary, the evaluations 
of K could be very complicated and its value very different, 
up to orders of magnitude toward much larger values. In this 
sense, we got the simplest and the lowest values for v — £ 
in this paper. 

Results shown throughout this paper were obtained 
working in the SPH framework. However, all results here 



discussed are not dependent on the adopted numerical code 
because physics of viscous dissipation is shown in its purely 
physical aspect, since we discussed a physical formulation 
for v. Parts of fluid, either as Lagrangian moving particles, 
or within Eulerian grid cells are conceptually the same both 
to the Euler and to the Navier-Stokes flow equations. In 
spite of fact that the artificial viscosity can also be tuned, 
according to the physical event considered, being SPH still 
an intrinsically artificially "viscous" technique, the evidence 
of successful results here shown demonstrates that, even in 
the worst conditions, v — £ correctly works without any 
ambiguity, free of any arbitrary parameters, also including 
those molecular characteristics that are uncommon in other 
formulations. 
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